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U.  S.  Naval  Postgraduate  School,  Monterey,  Calif. 

A  note  in  this  Journal,  Ref.  1,  discussed  the  problem  of  scheduling 
the  direction  p  of  constant  momentum  thrust  of  a  rocket,  which  loses 
mass  at  a  constant  rate,  so  that  it  transfers  to  a  known  earth  satel¬ 
lite  orbit  in  minimum  time  T  after  launching.  A  numerical  solution 
was  obtained,  using  rectangular  coordinates,  for  the  case  of  fixed 
launching  conditions.  The  method  of  Ref.  1  is  extended  here  to  solve 
the  problem  of  orbital  transfer  of  such  a  rocket  with  minimum  fuel 
consumption.  All  of  the  symbols,  units,  and  end  conditions  of  Ref.  1 
are  used  here  without  redefinition. 

Statement  of  the  Problem 

The  time  of  flight  T  in  minimum  fuel  transfer  must  be  longer  than 
in  the  minimum  time  transfer  of  Ref.  1,  unless  these  two  trajectories 
turn  out  to  be  identical.  This  implies  at  least  one  interruption  in 
rocket  thrust  during  minimum  fuel  transfer.  The  problem  solved  here 
assumes  exactly  one  such  interruption,  i.e.  launch  at  t=0,  thrust 
interruption  at  t=t^  ,  thrust  resumption  at  t=t2  ,  and  final  thrust 
termination  at  transfer  t=T.  The  problem  of  minimum  fuel  transfer  is 
equivalent  to  the  Lagrange  calculus  of  variations  problem  of  requir¬ 
ing  the  integral 

r  t 

J  =  JQ  (f+Xcp1+ucp2+ncp3+pcp4)dt  (1) 

to  be  stationary,  where  f  is  the  fuel  consumption  rate,  \,n,rr,p  are 
continua  of  Lagrangian  multipliers,  and  cp^  are  the  first 

order  equations  of  rocket  motion  of  Ref.  1.  The  f  function  and  the 


Received  September  ,  1963.  This  work  was  supported  by  the  Office 
of  Naval  Research.  Aid  of  the  Computer  Facility,  U.  S.  Naval  Post¬ 
graduate  School,  is  acknowledged. 

*  Professors,  Department  of  Mathematics  and  Mechanics. 

-1  - 


rocket  thrust' per  unit  remaining  mass  function  a  are  defined  as  fol¬ 
lows:  For  0<t<t.j  ,  f=1  and  a=cM/( 1 -Mt) g .  For  t^<t<t2  ,  f=0  and  a=0. 

For  t2<t<T,  f=1  and  a=cl^/[l  -ft(t+t.|  -^2)  ]g.  Note  that  for  t2<t<T 
3a/(3t^  =  -(3a/(3t2  =  ga  /c.  The  varied  time  subinterval  end  points  in 
Eq.(l)  are  taken  as  t^ +At^  ,  t2+At2  and  T+AT .  The  vanishing  first 
variation  6j  and  its  partial  integration  are  computed  as  in  Ref.  1. 
The  coefficients  of  <5u,6v,6x,6y  ,6p  in  6J=0  give  the  Euler  Eqs.(2) 
and  (3) ,  consisting  of  the  adjoint  equations 
\+tt=0  ,  U+p=0  , 

*+glxX+g2xM=0  ’  P+g1yx+g2y^=°  > 

and  the  control  equation 


(2) 


tan  p  =  pi/\  .  (3) 

The  coefficient  of  AT  in  6j=0  gives,  with  the  aid  of  Eq.(3),  the 
transversality  condition 

(a«A)T  =  (aA)T  =  1  (4) 

2  2l/2 

where  the  adjoint  vector  A  =  iX+ju,  A=|A|=(\+pi)'  ,  and 

a  =  a(i  cosp  +  j  sinp) .  The  coefficients  of  At^  and  At2  in  6J=0  give,  ^ 
with  the  aid  of  Eq.(4),  the  corner  conditions 

H(t1)  =  [aA]^  _  a2Adt  =  0,  H(t2)  =  0.  (5) 

Eqs.(5)  are  equivalent,  by  the  definition  of  a,  to 


A(t1 )  =  A(t2) 


and,  by  partial  integration,  to 

?T 

a^dt  =  0  . 


J XI  j 

s: 


(6) 

(7) 


Numerical  Solution 

Let  Xi jpu  , pi ,  i=1,2,3,4,  be  four  linearly  independent  solutions 
of  the  adjoint  Eqs.(2)  corresponding  to  the  columns  of  the  matrix  E(t) 
of  Ref.  1.  The  control  angle  p  of  Eq.(3)  is  defined  by 

tanp  =  (u1 +lpi2+m^3+n^4)/(^-| +1^2+m^3+n^4^  (8) 


2 


and  its  variation  6p  is  obtained  in  terms  of  6l,6m,6n  by  total  dif¬ 
ferentiation  as  in  Ref.  1 .  The  Bliss  fundamental  formulas  are  obtained 
by  assuming  that  a  solution  of  the  rocket  motion  equations  cp^ 
has  been  found,  corresponding  to  Eq.(8),  which  does  not  necessarily 
satisfy  the  terminal  conditions  at  t=T  or  the  corner  condition  Eqs.(5). 
Using  this  solution  and  holding  T  fixed,  but  allowing  t^  and  t2  to 

vary,  find  the  variation  of  the  vanishing  matrix  integral 

f  T 

JO  Cep-]  ,cp2  ,cp3  ,cp4]E( t)dt  =  0  (9) 

with  the  terminal  constraints  at  t=T  removed.  Since  the  columns  of 
E(t)  satisfy  the  adjoint  Eqs.(2),  one  obtains  the  system  of  Bliss 
formulas  in  the  1x4  matrix  equation 

[6u,6v,6x,6y]TE(T)  +  [G(t1 )-(apF)T]At1 

-  [G(t2)-(apF)T]At2  =  [0,6l,6m,6n]A  (10) 

where  the  matrix  A  has  been  defined  in  Ref.  1 ,  and  where  the  matrix 

G(t)  =  (apF) ^  a2pFdt  (11) 

where  the  2x4  matrix  F(t)  is  the  first  two  rows  of  E(t),  and  where 
the  matrix  p=[cosp,sinp] .  Substitution  of 

[6u,6v,5x,6y]T  =  [U-u,V-v,X-x,Y-y]T  +  [U-u,V-^,X-x,Y-y]TAT  (12) 
into  Eq.(lO)  gives  four  of  the  required  six  Newton-Raphson  equations 
for  the  determination  of  AT , At^ , At2 , 6l , 6m , 6n  on  the  varied  trajectory. 
The  remaining  two  equations  attempt  to  satisfy  the  corner  condition 
Eqs.(5)  on  the  varied  trajectory.  Involved  here  are  the  differentials 
da  =  6a  +  adt 

=  ( da/<3t.j )  At.j  +  (da/dt2)At2  +  (ga2/c)dt 
and  dA  =  6A  +  Adt 

=  [0 ,6l  jdmjdnjF'p*  +  (Xcosp+|i  sinp) dt 
where  the  primes  on  F  and  p  indicate  matrix  transposition.  Use  of 
Eqs.(6)  and  (14)  yields  the  Newton-Raphson  equation 

A(t1)At1  -;A(t2)At2  -  [0,6l,6m,6n][F*p*  ]t2  =  A(t2)  -  A(^)  .  (15) 


(13) 


(14) 
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Use  of  Eqs.(l3)  and  (14),  and  the  first  of  Eqs.(5),  yields  the 
Newton-Raphson  equation 

(aA)TAT  +  (K-aA)t  At1  -  K(t2)At2  +  [0 ,6l  ,6m,6n]G»  ( t1  )  =  -H(t.,)  (16) 


where 


K(t)  =  -  2(J)2^Ta3Adt 


(17) 


The  iteration  to  successive  varied  trajectories,  using  Eqs.(10), 
(12),  (15)  and  (16),  may  be  carried  out  as  in  Ref.  1.  Two  devices 
were  used  to  stabilize  the  course  of  the  iteration.  The  first  was  to 
adjust  the  m  and  n  values  of  the  new  T,t. ,t2,i,m,n  sextuple,  found  by 
solving  the  Newton-Raphson  equations,  to  satisfy  the  corner  condition 
Eqs.(5)  before  proceeding  with  the  next  iteration.  The  second  device 
was  to  modify  the  (_U,V,U,V,X,Y terms  in  the  Newton-Raphson  equations, 
before  solving  these  equations,  so  as  to  minimize  the  sum  of  the 
squares  of  the  elements  of  [U  -u ,  V-v  ,X-x ,  Y-yjrp. 

The  numerical  example  of  minimum  fuel  transfer  given  here  involves 
the  same  launching  conditions,  mass  loss  parameters  and  circular 
orbit  used  in  the  minimum  time  transfer  of  Ref.  1 .  The  results  for 
minimum  fuel  transfer  are  T=0. 353977,  t^=0. 210293,  t~=0. 275349, 
l=-0. 820196,  m=-0. 708727,  n=-1. 181390,  and  the  transfer  sector  angle 
B=0. 189345  rad.  Since  the  minimum  time  trajectory  of  Ref.  1  gave 
T=0. 289869,  the  net  fuel  saving  in  minimum  fuel  transfer  over  minimum 
time  transfer  is  measured  by  0. 289869-0. 353977+t0-t1  =  0.000948,  or 

Z  3 

+•  * 

an  unspectacular  one-third  per  cent.  Figure  1  shows  the  trajectories 
and  thrust  directions  for  minimum  time  and  minimum  fuel  transfer. 

The  semilogarithmic  plots  of  Fig.  2  show  the  different  behavior  of 
^ 'versus  time  in  the  two  problems.  For  some  reason  there  is  much  more 
difference  than  we  expect.  The  curve  increases  monOtonically  for  min¬ 
imum  time  transfer.  The  curve  for  minimum  fuel  shows  a  rather  char- 


i  I  1  .  >;•  ,1 
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acteristic  shape.  It  is  large  initially  and  decreasing;  then  it 
increases,  and  then  decreases.  If  the  final  decreasing  interval 
does  not  occur,  larger  values  of  T  lead  to  lower  values  of  fuel 
consumption,  as  may  be  partially  inferred  from  Eq.  (7). 

Reference 
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Bleick,  ¥.  E. ,  "Orbital  transfer  in  minimum  time,"  AIAA 
Journal  1,  1229-1231  (1963). 
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Fig.  1  Trajectories'  and  thrust  directions 
Fig.  2  A  versus  time 
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*  B  L  £  I C  K ( 3  3)  -  5  MINUTE  EXPRESS  7/31/63 

1 

PROGRAM  K IN FULL 

0 

CYVARSM)  =  XIJ  Y  VR  S  {  5 )  =XLM1  YVRS(9)  =XLM2  Y  VRS  (  13  )  =  XLM  3  Y  VRS  (  1  7  )  =XLM  4 

1 

C 

(25  =  X V  (6)  =  XMU  1  (10  5  =  XMU2 

(  1  4  )  =  XMU  3 

( 1 8) =XMU4 

2 

C 

(3)  -XX  (75  =XPI1  (  1 1  )  =XP I 2 

{ 1 5)=XPI 3 

(  1  9  )  =  X  P  I  4 

3 

c 

(4)  =  X  Y  (8)  =  X  R  0 1  (12)  =XR02 

{  1  6 ) =XRO  3 

(2C)=XRC4 

4 

CYVARS (21 ) = AA  1  2  YVRS { 25 ) = AA23  Y VRS ( 29 ) =AA44  YVRS ( 33 ) =C2A SG 

5 

C 

{ 2  2 )  =  A  A  1  3  (  26  5  -AA24  (30)=A2LAM 

(34) =C3ASG 

6 

c 

(23  5-AA14  (27 ! =  A  A  3  3  (31)=A3LAM 

(35  )=C4ASC 

7 

c 

(24)=AA22  (28)=AA34  (32)=C1ASC 

8 

2 

DIMENSION  YVARS (3 5 5 , AK(4,35 ) ,DY( 35),YC( 35) ,C (4) , XU (500) , 

XV ( 500  )  , 

9 

+  B( 6 ) , XX (  5  00 ) .XY (5  00 ) , TAU(  500)  , CAP  LAM { 500 ) 

tCGALAM (500) ,A2LAM(50C) , 

10 

+P1500) ,CC<4}, A{ 6, 6 ) ,AI (6, 6) ,DEL(6> ,CAPV(4) ,CAPVC(4) ,C IT1 (4 ), 

1  1 

+  CI T2 ( 4 ) ,CC  I  A2T2(4 )  , AAA( 4) , AA ( 4 ,4) , TVAR ( 6 ) 

, GCIA2T (4 ) 

12 

3 

ECU  I  VALENCE  (T , TV AR ( 1  )  ) » ( T 1 , T VAR ( 2 ) ) ,  ( T 2 , 

T V  AR ( 3 )  ), 

13 

1  (  EL, TVAR (4 ) >  ,  ( EM,TVAR( 5 ) ) ,  (  EN , 

TVAR (6  )  ) 

1  4  . 

4 

REARTB  =  20  925  OCC. 

1  5 

5 

GACCEL  =  32.086 

16 

6 

TUNIT  =  SQRTFtRE ARTH/GACCEL) 

17 

7 

CCC  -  10  000. 

18 

o 

V, 

COVERG  =  CCC/ (GACCEL  *  TUNIT) 

1  9 

9 

FMDOT  =  0.0036 

20 

10 

OMEGA  =  FMDOT  *  TUNIT 

21 

1  1 

VST  ART  =  0.585 

22 

12 

THETA  =  0.928 

23 

13 

T  -  0.353  966  649 

24 

14 

T 1  =  C.210  274  52 C 

25 

1  5 

T 2  -  0.275  321  127 

26 

16 

EL  =-0.320  214  924 

27 

17 

EM  =-0.708  762  302 

28 

18 

EN  =-1.181  45o  519 

29 

19 

R  =  1.075  693  925 

30 

20 

V  =  S  G  R  T  F ( 1.0/R) 

31 

21 

VSQOR  =  V*V/R 

32 

22 

CB  =  0.189  335  935 

33 

23 

TFIN  =  0.28972  53036 

34 

24 

XSTEP  =  TFIN/ 116. 

35 

25 

XU(1)  =  V ST ART  *  CCSF(THETA) 

36 

26 

XV (1)  =  VST  ART  *  SINE (THETA) 

37 

27 

XX ( 1  !  =  C  .0 

38 

28 

XY ( 1 )  =  1.0 

39 

29 

T  AU  (  1  )  '=  0.0 

40 

31 

A2L AM (15=  O.C 

42 

32 

C  (  1  )  =  C  .  0 

43 

33 

C ( 2  )  =  C.5 

44 

34 

C ( 3 )  =  0.5 

45 

35 

C  (  4  5  =  1.0 

46' 

KK  =  C 

46. 1 

36 

DO  27  1  L= 1 s 3 

47 

57 

XV A R  =  0.0 

48 

38 

YVA'RS  (  1  )  =  XU  (  1  ) 

49 

39 

YVARS ( 2 )  =  XV ( 1 ) 

50 

40 

YVARS ( 3  )  =  XX ( 1  ) 

51 

4  1 

YVARS (4)  =  X  Y {  1  ) 

52 

42 

C  A  P  L  A  M (  1  ) =  SORT F(  1.0  +  EL*EL ) 

53 

425 

P {  1  )  =  57.2957  *  ATANF(EL) 

53.5 

45 

XA  =  COVERG  *  OMEGA 

54 

44 

CGALAH 1 )  =  COVERG  *  X A  *  CAPLAM(l) 

55 

45 

DO  4  6  I  =6  »  35 

56 

46 

YVARS (I)  =  0.0 

57 

47 

DO  48  1  =  5,  20,5 

58 

48 

YVARSUi  =  1.0 

59 

49 

N 1  =  11 4X STEP  +  1  .0 

60 

50 

XN1  =  N 1 

61 

51 

STEP!  =  T1/XN1 

62 

52 

N2  =  (T2-T 1  ) /XSTEP  +  1.0 

63 

53 

X  N 2  =  N2 

64 

54 

STEP2  =  ( T2-T  1  )  /XN2 

* 

65 

55 

N2  =  M  +  N2 

66 

56 

N5  =  ( T- T  2 ) /X  STEP  +  1.0 

67 

57 

XN3  =  N 3 

68 

58 

STEPS  =  ( T-T2 ) /XN  3 

69 

59 

N3  =  N2  +  N3 

70 

60 

SINP  =  S  IMF ( Be  > 

7  1 

61 

CO  SB  =  CU  Sr { D  B ) 

72 

62 

CAPV (  1  )  =  V  *  COS  6 

73 

63 

CAPV(  25  =  -V  *  SI  N 6 

74 

64 

CAPV ( 3)  =  R  *  S IN  B 

75 

65 

C A P V  (  4  )  =  R  *  COS  B 

76 

66 

CAPVD(l)  =  -VSQDR  *  SINB 

77 

67 

C AP VD ( 2  5  =  -VSGDR  *  COSB 

78 

-  6  - 


68 

CAPVD ( 3  )  =  CAPV ( 1  ) 

79 

69 

CAPVDI u )  =  CAPV (2 ) 

80 

70 

N4  =  N3  +  1 

81 

71 

CO  206  K=2,N4 

82 

72 

IF  (  N  1  +  1  -  K  5  75,73,73 

83 

73 

STEP  =  STEP1 

84 

74 

GC  T 3  79 

85 

75 

IF  (N2+1-K)  78,76,76 

86 

76 

STEP  =  STEP2 

87 

77 

GO  TC  79 

88 

76 

STEP  =  STEP3 

89 

79 

DC  124  1=1,4 

90 

80 

XC  =  XVAR  +  C ( I )  *  STEP 

9  1 

81 

DC  8  2  J  =  1 ,35 

92 

82 

YC ( J  )  =  YVARS(J)  +  C (  I )  *  AK  (  I  —  1  ,  J  ) 

93 

83 

XL  AM  =  YC ( 5  )  +  EL  *  Y  C ( 9 )  +  £M«YC(13)  +  EN*YC(17) 

94 

eu 

XMU  =  YC ( 6 )  +  EL  * YC (  10)  +  EM*YC(14)  +  EN*YC(18) 

95 

85 

CLAM  =  SCRTF ( XL AM* *2  +  XMU**2) 

96 

86 

COSP  =  X  L  AM/CLAM 

97 

87 

SIMP  =  XMU /CL  AM 

98 

88 

IF  (Nl+l-K)  91,89,89 

99 

89 

XA  =  C  0  V  L  R  G  *  CMEGA/d.O  -  CMEGA  *  XC  ) 

100 

90 

GO  TO  95 

101 

91 

IF  (N2+ 1-K )  94 ,92,92 

102 

92 

XA  =  C.C 

103 

93 

GO  TO  95 

1  04 

9  4 

XA  =  CCV  2RG*0M EGA / (  1  .  C- OMEGA -»  { XC-T2  +  T 1 )  ) 

105 

95 

XR  =  SCRTF(YC(3)**2  +  YC(4)**2) 

106 

96 

CY ( 1 )  =  -YC(3  )/XR**3' +  XA*COSP 

107 

97 

DY ( 2 )  =  -YC(4)/XR**3  +  XA*SINP 

108 

98 

DY ( 3  )  =  YC (  1 ) 

1  09 

99 

DY  (  4  )  =  YC  (  2  ) 

1  10 

ICC 

G 1 X  =  (2.*YC(3)**2  -  YC(4)**2)/XR**5 

1  1  1 

101 

G1Y  =  3.«YC(3)*YC (4)/XR**5 

1  12 

102 

G2X  =  G1Y 

1  13 

1  03 

G2Y  =  (2.*YC(4)»*2  -  YC (3 )**2)/XR**5 

1  1  4 

104 

DO  108  M  =5 ,17,4 

1  15 

105 

DY  (  M  )  =  -  Y  C  ( M  +  2  ) 

1  16 

106 

0Y( 1  )  =  -YC ( M+3  ) 

1  17 

107 

DY  (  M  +  2  )  =  - G  1  X <•  YC  (  M  )  -  G2X*YC(M+1) 

1  18 

108 

DY  (  M  +  3  )  =  -  G  1  Y  *  YC  (  M  )  -  G2Y*YC(M  +  1) 

1  19 

109 

DO  11C  M  =  1  ,  4 

120 

1  10 

A A A ( M  )  =  COSP*YC(  4*M+2)  -  SI NP* YC ( 4*M  +  1  ) 

121 

1  1  1 

DO  11 2  M=  1 , 3 

122 

112 

DY ( 2 0  +  M  )  =  XA*AAA ( 1  ) *AAA( M+l ) /CLAM 

123 

1  13 

DO  114  M  = 1  ,3 

124 

1  14 

C  Y ( 2  3  +  M )  =  XA*AAA  (2)*AAA(M+1  l/CLAM 

125 

1  15 

DY ( 27 )  =  XA*AAA (3 ) *AAA( 3 ) /CLAM 

126 

1  16 

DY ( 2  8  )  =  XA*AAA(3  )*AAA( 4) /CLAM 

1  27 

1  17 

DY ( 29  )  =  XA*AAA(4  )*AAA(4)/CLAM 

123 

1  18 

DY  (-30  )  =  XA*XA*CL  AM 

129 

1  19 

DY  (  3  1  )  =  XA»DY(  30  ) 

1  30 

120 

DO  122  M=l,4 

131 

121 

CC  (  M  )  =  YC  (  4*M+1  )  *COSP  +  YC(4*-M+2)*SINP 

1  32 

122 

DY (31+  K )  =  CC(M)*XA*XA 

133 

123 

DO  124  J  =  1  ,  3 5 

1  34 

124 

AK(  I  ,  J  )  =  STE P*-DY  (  J  ) 

135 

125 

DO  126  J=  1 ,35 

136 

126 

127 

128 

129 

1 30 

131 

132 

133 

134 

135 
1  36 
137 
1  38 
1  39 
140 
1  46 
1  47 

148 

149 
1  50 

151 

152 

153 


YVARS ( J )  =  YVARS( J ) 

XV  AR  =  XVAR  +  STEP 
TAU(K)  =  TAU(K-l)  + 

XU  (  K  )  =  WARS  (  1 ) 

XV  (K)  =  WARS  (2) 

XX ( K  )  =  YV ARS ( 5 ) 

XY ( K )  =  Y VARS ( 4 ) 

XL  AM  =  YVARS ( 5  )  + 

XMU  =  YVARS ( 6 )  + 

CLAM  =  SCRTF { XLAM**-2 
CAPLAN(K)  =  CLAM 
COSP  =  XL AM/CLAM 
S1NP  =  X MU/CLAM 
DC  1  4  C  M  =  1 , 4 
CC ( M )  =  YVARS t 4*M+ 1 )*C0SP 
DO  148  M  =  5  »  1  7 , 4 
D  Y  (  M  )  =  -  WARS  (  M+  2  ) 

DY { M+ 1 )  =  -YVARSt  M+3 ) 
CGALAM(K)  =  CCV  ER  G*X  A*CL  AM 
A2 L  AM  (  K  )  =  YVARS(  30 
IF  ( X L AM  )  1  59,152,159 

IF  (XMU)  1  57,  15  5,1  53 
P(K)  =  90. 0 


+  (  AK(1,J)+2.*AK(2,J)+2.*AK(3,J)+AK(4,J)  )/6. 

STEP 


E  L  *Y VAR S ( 9  )  +  EM*YVARS(13)  +  EN*YVARS(17) 
EL*  YVAR  S (10)  +  E  M* YVARS (14)  +  EN*YVARS(18) 
+  XMU* *2 ) 


+  YVARS(4*M+2)*SINP 


-  7  - 


1  3  7 
138 
1  39 

140 

141 

142 

143 
1  44 

145 

146 

147 

148 

149 
1  50 

15  1 

157 

158 

159 

160 

16  1 
1  62 

163 

164 


154  GO  TC  165 

155  P(K)  =  0.0 

156  GO  TO  165 

157  P(K)  =  -00.0 

158  GO  TO  165 

159  P ( K )  =  57.2957  *  A  T  ANF ( XMU/XLAM) 

160  IF  (XLA.M)  161,165,165 

161  IF  C  XT-'U  >  164,  16  2,lo2 

162  P  (  K )  =  P ( K )  +  180 .0 

163  GO  TO  165 

164  p ( K )  =  P ( K )  -  180. C 

165  X  L  A  M OCT  =  D Y ( 5 )  +  EL  *DY I 9 )  +  EM*DY(13)  +  EN»DY(17> 

166  XMUOOT  *  D Y ( 6 )  +  EL*DY(10)  +  EM*DY(14)  +  EN*DY{18) 

167  IF  (Nl+l-K )  1  85,  168,  183 

168  00  169  M=  1 , 4 

169  CIT1 IN)  =  CC( M) 

170  A2LAMT1  =  XA*XA*OLAM 

171  CLANDT1  =  COSP*XLAKCCT  +  SINP*XMUDOT 

172  CGAL0T1  =  CCVERG*-XA*CLAMDT1 

173  CLAMT1  =  CLAM 

174  C2ACGT1  =  CCVERG*XA*CC(2) 

175  C3ACGT1  =  COVER G*XA»CC( 3) 

176  C4ACGT1  =  CCVERG*  XA*CC I  4 ) 

177  ATI  =  XA 

178  CGALMT1  =  COVERG*  XA*CLAM 

179  0A2LMT2  =  YVARS(3C) 

1  SC  CA3LMT2  =  2.0*YVARS(31 J/COVERG 

181  DO  182  M  =  1  ,4 

182  QC I A2T2 ( M )  =  YVARSI31+M) 

183  IF  (M2+1-K)  190, 184, 190 

184  00  185  N  =  1  , 4 

185  C  I T2  I  M  )  =  CC(M) 

186  A2LAMT2  =  A  T 1  *  A  T 1  *  CLAM 

187  CGALAM(K)  =  COVER G*A T  1  *C LAM 

188  CLAMDT2  =  CCSP*XLAMDOT  +  SINP*XMUDOT 

189  CLAPT2  =  CLAM 

190  IF  ( N  4  —  K )  206,  19  1,206 

191  CLAMOT  =  COSP^XLA  MOOT  +  SINP*XMUDOT 

192  CGALOT  =  COVERG*X  A»CLAMDT 

193  C2ACGT  =  COVERG*X A*CC ( 2) 

194  C3ACGT  =  COVERG*XA*CC ( 3) 

195  C4ACGT  =  COVERG*X A*CC (4) 

196  CGALMT  =  COVERG <XA*CLAM 

197  QA2LMT  =  YVARSI 30 ) 

198  CA3LMT  =  2 . 0* YVAR S ( 3 1 ) /COVERG 

199  DO  20C  M  = 1 , 4 

20C  QCIA2T(M)  =  Y VARS ( 3 1+M ) 

201  XR  =  SQRTF( YVARSI  3)**2  +  YVARS(4)**2) 

202  CY  (  1  )  =  -YVARSI  3)  /XR*-*3  +  XA*COSP 

203  DY I  2  )  =  -YVARSI 4)  /XRt*3  +  XA*SINP 

204  DY  ('3  )  =  YVARSI  1  ) 

205  DY ( 4 )  =  YVARSI 2 ) 

206  CONTINUE 

207  PRINT  203 

2C8  FORMAT! 1 PG6X2HEL1 3X2HEM1 3X2HEN1 3X2HBB 1 3X2HT 11 3X2HT21 3X1 HT ) 

209  PRINT  2 1 C ,  EL , EM, EN , BB, T 1 , T2 , T 

210  FORMAT  (7F15.9) 

211  PRINT  212 

212  FORMAT  (  1 HC6X 2HN1  1 3X 2 HN2 1 3X2 HN3 1 3X 1  HU  1 4 X 1 H V 1 4 X 1 HXl  4X 1 H Y ) 

213  PRINT  214,  N1,N2, N3, XUI N4 ) , XV (N4 ) , XX  I N4 ) , X Y I N4 ) 

214  FORMAT  I  3  I  1 5 , 4 F 1 5 . 7 ) 

215  PRINT  216 

216  FORMAT  I  1H05X4HCAPU1 1 X4HCAPV1 1 X4HCAPX1  1  X4HCAPY) 

217  PRINT  218,  (CAPVIM),  M  =  1 , 4 ) 

218  FORMAT  (4F15.7) 

219  DC  224  1=1,4 

220  A  I  I , 1  )  =  0.0 

221  BII)  =  C.O 

222  DO  224  J=l,4 

223  A  I  I ,  1  )  =  A I  1 ,  1  )  +  YVARS(4*I+J )*(DY( JJ-CAPVDI J) ) 

224  BII)  =  BII)  +  YVARS(4*I  +  J)MCAPV(J)-YVARS(J)  ) 

225  DO  227  1=1,4 

226  All, 2)  =  AT1*CIT1(I)  +  ( OC I A2 T I  I ) - CC I A2 T2 I  I ) ) /COVE RG 

227  All, 3)  =-AT  1*C I T2 I  I )  -  I QCI A2T I  I ) -QC I A2T2 I  I )  ) /COVE RG 

228  DO  23C  J=2,4 

229  A  A I  1  , J ) =  YVARSI 19+J ) 

2 30  AA I  2 , J ) =  YVARS I  22  + J ) 

231  AA I  3, 3 ) =  YVARS I  27 ) 

232  A  A ( 3 , 4 ) =  YVARSI28 ) 

233  A A ( 4 , 4 ) =  YVARSI 29 ) 

234  DO  23c  1=1,4 


1  65 
166 

167 

168 
1  69 

170 

171 

172 

173 

174 

175 

176 

177 

178 
1  79 
1  80 
181 
182 
1  33 

184 

185 

186 
1  87 
188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 
207 
20  8 

209 

210 
21  1 
212 

213 

214 
21  5 
216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 

244 

245 
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245  A  (  5 , 6  ) 

246  B  (  5  ) 

247  A( 6,  1  ) 

248  A  (  6 , 2  ) 


235  DO  236  J  =  1  ,  I 

236  AA( I, J)  =  A  A  (  J  , I) 

237  CO  239  1  =  1  ,4 

238  DO  239  J=1 , 3 

239  A  (  I  ,  J  +  3  )  =  A  A  (  I,J  +  1) 

240  A( 5, 1  )  =  -CGALDT 

241  A  (  5  >  2  )  =  A2LAMT1  +  CGALDT 1  +  C  A  3  L  M  T  -  CA3LMT2  -  XA*XA*CLAM 

242  A  (  5 , 3  )  =-A2LAMT2  -( CA3LMT  -  CA3LMT2)  +  XA*XA*CLAM 

243  A ( 5  » 4  )  =  QC  I A2T ( 2  )  -  CCI42T2(2)  +  C2ACGT1  -  C2ACGT 

244  A  (  5 1 5 )  =  QCIA2K3)  -  QCIA2T2(3)  +  C3ACGT1  -  C3ACGT 

“  “  CCIA2K4)  -  QC  I  A  2T2  (  4  )  +  C4ACCT1  -  C4ACGT 

CGALMT  -  CGALMT1  -  QA2LMT  +  QA2LMT2 
0.0 

CLAMDT1 

249  A ( 6 , 3  )  =-CL AMDT2 

250  DO  251  J  =  2,4 

251  A ( 6  »  J  +2  )  =  CITl(J)  -  C  I T2 ( J ) 

252  H ( 6 )  =  CLAMT2  -  C  L AMT  1 

272  DO  274  I  =  1  ,N4 

273  CGALAM(I)  =  CGALA  M ( M4 )  -  CGALAM(I) 

274  A 2 L A M (  I  )  =  A2LAM{ N4 )  -  A2 LAM (  I  ) 

IF  (KK)  320,320,300 

30C  KK  =  KK- 1 

DET  =  A(5,5)«A(6, 6)-A(5,6)*A(6,5) 

EM  =  EM  +  ( B( 5)*A (6,6)-B(6)*A(5,6) J/DET 
EN  =  EN  +  (R(6)*A(5,5)-B(5)*A(6,5)  )/DET 
PRINT  301 

301  FORMAT  (  1 H04X3HTA U8X6HCAPLAM7X6HCGALAM8X5HA2L AM/ ) 

PRINT  302  ,  {  T  AU  (  I  )  ,C’APLAM(  I  )  ,CGALAM(  I  )  ,A2LAM(  I  )  ,  1  =  1, N4) 

302  FORMAT  (4F13.9) 

GO  TO  37 

320  Dbl  =  ( XU( N4 )-CAP V( 1 ) )/CAPV( 2) 

DB2  =-( XV ( N4) -CAP V( 2 ) ) /CAPV ( 1 ) 

DB3  =  ( XX ( N 4 ) -C  AP  V ( 3 )  ) /C  AP  V ( 4 ) 

DB4  =-(XY(N4)-CAPV(4) )/CAPV( 3) 

BB  =  BB  +  ( DB1+DB2+DB3+DB4) /4. 

SINB  =  SINF(BB) 

CCS  D  =  CCSF(BB) 

CAPV( 1 )  =  V*COSB 
CAP V( 2 )  =-V*S  INB 
C A  PV (  3 )  =  R»S INB 
CAPV ( 4 )  =  R  *C  CSB 
CAP  VD ( 1 ) =- VSQDR*S INB 
CAP VD ( 2 ) =- V  SGDR*C  OS  B 
CA  PVD ( 3 ) =  C  AP  V (  1) 

CAPVD ( 4 ) =  C A  P V ( 2) 

PRINT  208 

PRINT  2 1 C ,  EL, EM, EN, BB,T1  ,T2,T 

321  PRINT  216 

322  PRINT  218,  (CAPV(M),  M=l,4) 

323  DO '326  1=1,4 

324  A ( I , 1 )  =  0.0 

325  B (  I )  =  0.  C 

326  DO  32  8  J  =  1  ,  4 

327  A (  I ,  1  )  =  A ( I ,  1  )  +  YVARS(4*I  +  J)*(DY( J)-CAPVD(  J)  ) 

328  B(I)  =  B (  I  )  +  YVARS(4*I+J)*(CAPV(J)-YVARS(J)  ) 


253  CALL  GAUSS3  (6,  0.1E-09,  A,  AI ,  KER) 

254  PRINT  255,  KER 

255  FORMAT  ( 5H0KER=I1  ) 

256  IF  (KER-2)  258,257,258 

257  STOP  257 

258  DO  261  1=1,6 

259  DEL (I)  =  0.0 

260  CO  261  J=  1  »  6 

261  DEL (  I  )  =  DEL (  I  )  +  AI(I,J)*B(J) 

262  DO  263  1=1,6 

263  TVAR(I)  =  TVARII)  +  DELI  I  ) 

264  BB  =  BB  +  V*DEL ( 1 ) /R 

265  PRINT  208 

266  PRINT  210,  EL  ,  EM,  EN  ,  RB  ,  T  1  ,  T2  ,  T 

267  IF  (T2-T)  303,303,307 

303  IF  ( T 1  —  T )  304,304,307 

304  IF  I  T  2  )  307,305,  305 

305  IF  ( T 1 )  307, 306, 306 

306  IF  (T1-T2)  271,271,307 

307  GO  TO  275 
271  CONTINUE 

275  PRINT  276 

276  FORMAT!  1 H04X3HTAU 1 0X2HXU 1 1X2HXV1 1 X2HXX 1 1 X2HXY9X6HC 
1X5HA2LAM9X1HP/) 

278  PRINT  279,  ( TAUI I  )  ,  XU ( I ) , XVI  I  ) , XX ( I )  , XY ( I ) , CAPLAMI 

-  9  - 


246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 
263 

283 

284 

285 
263.  1 

263.2 

263.3 

263. 4 
263.  5 

263.6 

263.7 
263.90 
263.92 
263.95 


263.99 

263991 


264 

26  5 
266 

267 

268 

269 

270 

27  1 

272 

273 

274 
•  275 

276 

277 
2  78 
278.  1 

279 
279.1 

280 
280.  r 
282 
286 

APLAM7X6HCGALAM8287 

288 

I) ,CGALAM( I ) ,  289 


i 


A2LAM  I  )  ,  P  (  I  )  ,  1=1 , N4) 

FORMAT  (SF1 3.9,  F  13.2  ) 

290 

;79 

291 

!BC 

STOP  280 

292 

>81 

2ND 

293 

SURROLT I NE  GAUSS3  (N, E P , A , X , K E R ) 

294 

DIMENSION  A ( 6 , 6 ) ,  X { 6 , 6 ) 

295 

CO  1  1  =  1  ,  N 

C0030 

CO  1  J=1,N 

0004  C 

1 

X(  I, J  )  =  C.C 

000  5C 

DC  2  K=1,N 

00C6C 

c 

X  (  K  ,  K  )  =  1  .C 

00070 

10 

DO  34  L  =  1  ,  N 

ocoec 

KP  =  C 

000  9  C 

Z  =  0.0 

OC  100 

DO  12  K  =  L »  N 

00110 

IF(Z-ARSF(A(K,L)) ) 11 , 1 2, 1 2 

00  120 

M 

Z=APSF { A{ K, L) ) 

001  3C 

KP  =  K 

00140 

12 

CONTINUE 

00  1  50 

IF ( L-KP ) 13,20,20 

00160 

13 

CC  14  J=  L  ,  N 

OC  170 

Z  =  A (L, J  ) 

ooiao 

A(L,J )  =  A { K  P , J ) 

00190 

lit 

A{ KP, J )=Z 

0 02  OC 

CO  15  J= 1 , N 

002  10 

Z  =  X ( L , J  ) 

00220 

X(L, J ) =X ( KP , J ) 

0023C 

15 

X  (  K  P  ,  J  )  =  Z 

0024C 

20 

IFIABSFI A(L,L) )-EP) 50,50,30 

00250 

iO 

I F  (  L— N ) 3 1 ,34,34 

00260 

51 

LP  1  =  L  +  1 

C0270 

CC  36  K  =  LP  1  »  N 

OC  230 

I  F  (  A(  K,  L  )  )  32, 36,32 

00290 

32 

RATIO=A{K,L)/A(L,  L) 

0C300 

DO  33  J=LP1,N 

0031  C 

33 

A(K,J)=A(K,J ) -RAT IO*A(L, J) 

0C320 

DO  35  J=1 ,N 

00330 

35 

X(K,J  )=X(K,J)-RATIC*X(L,J) 

00340 

36 

CONTINUE 

00350 

34 

CONTINUE 

00360 

40 

DO  43  1=1, N 

C0370 

I  I  =N+  1- I 

0038C 

CO  43  J=1 , N 

0C390 

S=0 .0 

0040C 

I F (  I  I - N ) 4 1 , 43,43 

0041  C 

41 

I  I  P  1  =  1 1  +  1 

0  042  0 

DO  42  K= I  IP  1  ,  N 

00430 

42 

S=  S  +  A (  I I ,K )*X (K,J  ) 

00440 

43 

X{  I  I , J )  =  { X (  II,J)— S)/A(  II  ,11) 

00450 

K  ER=  1 

00460 

50 

RETURN 

00470 

KER  =  2 

00480 

END 

0049C 

END 

296 
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